home *** CD-ROM | disk | FTP | other *** search
/ 3D GFX / 3D GFX.iso / amiutils / i_l / irit5 / triv_lib / geomat4d.c next >
C/C++ Source or Header  |  1995-12-30  |  8KB  |  201 lines
  1. /******************************************************************************
  2. * GeoMat4d.c - Trans. Matrices , Vector computation, and Comp.geom, in 4D.    *
  3. *******************************************************************************
  4. * Written by Gershon Elber, November 1994.                                    *
  5. ******************************************************************************/
  6.  
  7. #include <math.h>
  8. #include <stdio.h>
  9. #include "irit_sm.h"
  10. #include "triv_loc.h"
  11.  
  12. /*****************************************************************************
  13. * DESCRIPTION:                                                               M
  14. * Computes a hyperplane in four space through the given four points.         M
  15. * Based on a direct solution in Maple of:                                    M
  16. *                                         M
  17. *  with(linalg);                                 V
  18. *  readlib(C);                                     V
  19. *                                           V
  20. *  d := det( matrix( [ [x - x1, y - y1, z - z1, w - w1],             V
  21. *                      [x2 - x1, y2 - y1, z2 - z1, w2 - w1],               V
  22. *                      [x3 - x2, y3 - y2, z3 - z2, w3 - w2],             V
  23. *                      [x4 - x3, y4 - y3, z4 - z3, w4 - w3]] ) );         V
  24. *  coeff( d, x );                                 V
  25. *  coeff( d, y );                                 V
  26. *  coeff( d, z );                                 V
  27. *  coeff( d, w );                                 V
  28. *                                         *
  29. * PARAMETERS:                                                                M
  30. *   Pt1, Pt2, Pt3, Pt4:  The four points the plane should go through.        M
  31. *   Plane:               Where the result should be placed.             M
  32. * RETURN VALUE:                                                              M
  33. *   int:     TRUE if successful, FALSE otherwise.                            M
  34. *                                                                            *
  35. * KEYWORDS:                                                                  M
  36. *   TrivPlaneFrom4Points, hyper plane, plane                                 M
  37. *****************************************************************************/
  38. int TrivPlaneFrom4Points(TrivPType Pt1,
  39.              TrivPType Pt2,
  40.              TrivPType Pt3,
  41.              TrivPType Pt4,
  42.              TrivPlaneType Plane)
  43. {
  44.     Plane[0] = - Pt2[1] * Pt3[3] * Pt4[2] -
  45.                  Pt1[1] * Pt2[2] * Pt3[3] +
  46.                  Pt2[1] * Pt3[2] * Pt4[3] -
  47.                  Pt1[1] * Pt3[2] * Pt4[3] +
  48.                  Pt1[1] * Pt2[2] * Pt4[3] +
  49.                  Pt1[1] * Pt3[3] * Pt4[2] -
  50.                  Pt1[1] * Pt2[3] * Pt4[2] +
  51.                  Pt1[1] * Pt2[3] * Pt3[2] -
  52.                  Pt3[1] * Pt2[2] * Pt4[3] +
  53.                  Pt3[1] * Pt1[2] * Pt4[3] +
  54.                  Pt3[1] * Pt2[3] * Pt4[2] -
  55.                  Pt3[1] * Pt1[3] * Pt4[2] -
  56.                  Pt2[1] * Pt1[2] * Pt4[3] +
  57.                  Pt2[1] * Pt1[2] * Pt3[3] +
  58.                  Pt2[1] * Pt1[3] * Pt4[2] -
  59.                  Pt2[1] * Pt1[3] * Pt3[2] +
  60.                  Pt4[1] * Pt2[2] * Pt3[3] -
  61.                  Pt4[1] * Pt1[2] * Pt3[3] +
  62.                  Pt4[1] * Pt1[2] * Pt2[3] -
  63.                  Pt4[1] * Pt2[3] * Pt3[2] +
  64.                  Pt4[1] * Pt1[3] * Pt3[2] -
  65.                  Pt4[1] * Pt1[3] * Pt2[2] -
  66.                  Pt3[1] * Pt1[2] * Pt2[3] +
  67.                  Pt3[1] * Pt1[3] * Pt2[2];
  68.     Plane[1] = - Pt2[0] * Pt3[2] * Pt4[3] +
  69.                  Pt2[0] * Pt3[3] * Pt4[2] -
  70.                  Pt1[0] * Pt2[2] * Pt4[3] +
  71.                  Pt1[0] * Pt3[2] * Pt4[3] +
  72.                  Pt1[0] * Pt2[2] * Pt3[3] -
  73.                  Pt1[0] * Pt3[3] * Pt4[2] +
  74.                  Pt1[0] * Pt2[3] * Pt4[2] -
  75.                  Pt1[0] * Pt2[3] * Pt3[2] +
  76.                  Pt3[0] * Pt2[2] * Pt4[3] -
  77.                  Pt3[0] * Pt2[3] * Pt4[2] -
  78.                  Pt3[0] * Pt1[2] * Pt4[3] +
  79.                  Pt3[0] * Pt1[3] * Pt4[2] +
  80.                  Pt2[0] * Pt1[2] * Pt4[3] -
  81.                  Pt2[0] * Pt1[2] * Pt3[3] -
  82.                  Pt2[0] * Pt1[3] * Pt4[2] +
  83.                  Pt2[0] * Pt1[3] * Pt3[2] -
  84.                  Pt4[0] * Pt2[2] * Pt3[3] +
  85.                  Pt4[0] * Pt2[3] * Pt3[2] +
  86.                  Pt4[0] * Pt1[2] * Pt3[3] -
  87.                  Pt4[0] * Pt1[3] * Pt3[2] -
  88.                  Pt4[0] * Pt1[2] * Pt2[3] +
  89.                  Pt4[0] * Pt1[3] * Pt2[2] +
  90.                  Pt3[0] * Pt1[2] * Pt2[3] -
  91.                  Pt3[0] * Pt1[3] * Pt2[2];
  92.     Plane[2] =   Pt2[0] * Pt3[1] * Pt4[3] -
  93.                  Pt2[0] * Pt4[1] * Pt3[3] -
  94.                  Pt1[0] * Pt2[1] * Pt3[3] -
  95.                  Pt1[0] * Pt3[1] * Pt4[3] +
  96.                  Pt1[0] * Pt2[1] * Pt4[3] +
  97.                  Pt1[0] * Pt4[1] * Pt3[3] -
  98.                  Pt1[0] * Pt4[1] * Pt2[3] +
  99.                  Pt1[0] * Pt3[1] * Pt2[3] -
  100.                  Pt3[0] * Pt2[1] * Pt4[3] +
  101.                  Pt3[0] * Pt4[1] * Pt2[3] +
  102.                  Pt3[0] * Pt1[1] * Pt4[3] -
  103.                  Pt3[0] * Pt4[1] * Pt1[3] -
  104.                  Pt2[0] * Pt1[1] * Pt4[3] +
  105.                  Pt2[0] * Pt1[1] * Pt3[3] +
  106.                  Pt2[0] * Pt4[1] * Pt1[3] -
  107.                  Pt2[0] * Pt3[1] * Pt1[3] +
  108.                  Pt4[0] * Pt2[1] * Pt3[3] -
  109.                  Pt4[0] * Pt3[1] * Pt2[3] -
  110.                  Pt4[0] * Pt1[1] * Pt3[3] +
  111.                  Pt4[0] * Pt3[1] * Pt1[3] +
  112.                  Pt4[0] * Pt1[1] * Pt2[3] -
  113.                  Pt4[0] * Pt2[1] * Pt1[3] -
  114.                  Pt3[0] * Pt1[1] * Pt2[3] +
  115.                  Pt3[0] * Pt2[1] * Pt1[3];
  116.     Plane[3] = - Pt2[0] * Pt3[1] * Pt4[2] +
  117.                  Pt2[0] * Pt4[1] * Pt3[2] +
  118.                  Pt1[0] * Pt2[1] * Pt3[2] +
  119.                  Pt1[0] * Pt3[1] * Pt4[2] -
  120.                  Pt1[0] * Pt2[1] * Pt4[2] -
  121.                  Pt1[0] * Pt4[1] * Pt3[2] +
  122.                  Pt1[0] * Pt4[1] * Pt2[2] -
  123.                  Pt1[0] * Pt3[1] * Pt2[2] +
  124.                  Pt3[0] * Pt2[1] * Pt4[2] -
  125.                  Pt3[0] * Pt4[1] * Pt2[2] -
  126.                  Pt3[0] * Pt1[1] * Pt4[2] +
  127.                  Pt3[0] * Pt4[1] * Pt1[2] +
  128.                  Pt2[0] * Pt1[1] * Pt4[2] -
  129.                  Pt2[0] * Pt1[1] * Pt3[2] -
  130.                  Pt2[0] * Pt4[1] * Pt1[2] +
  131.                  Pt2[0] * Pt3[1] * Pt1[2] -
  132.                  Pt4[0] * Pt2[1] * Pt3[2] +
  133.                  Pt4[0] * Pt3[1] * Pt2[2] +
  134.                  Pt4[0] * Pt1[1] * Pt3[2] -
  135.                  Pt4[0] * Pt3[1] * Pt1[2] -
  136.                  Pt4[0] * Pt1[1] * Pt2[2] +
  137.                  Pt4[0] * Pt2[1] * Pt1[2] +
  138.                  Pt3[0] * Pt1[1] * Pt2[2] -
  139.                  Pt3[0] * Pt2[1] * Pt1[2];
  140.  
  141.     Plane[4] = -(Plane[0] * Pt1[0] + Plane[1] * Pt1[1] +
  142.          Plane[2] * Pt1[2] + Plane[3] * Pt1[3]);
  143.  
  144.     return ABS(Plane[0]) > EPSILON ||
  145.        ABS(Plane[1]) > EPSILON ||
  146.        ABS(Plane[2]) > EPSILON ||
  147.        ABS(Plane[3]) > EPSILON;
  148. }
  149.  
  150. #ifdef TEST_PLANE_4D
  151. void main(void)
  152. {
  153.     int i, j;
  154.     TrivPlaneType Plane;
  155.     TrivPType
  156.     Pt1 = { 1, 0, 0, 0 },
  157.     Pt2 = { 0, 1, 0, 0 },
  158.     Pt3 = { 0, 0, 1, 0 },
  159.     Pt4 = { 0, 0, 0, 1 };
  160.  
  161.     TrivPlaneFrom4Points(Pt1, Pt2, Pt3, Pt4, Plane);
  162.     printf("[%lf %lf %lf %lf %lf]\n",
  163.        Plane[0], Plane[1], Plane[2], Plane[3], Plane[4]);
  164.  
  165.     TrivPlaneFrom4Points(Pt1, Pt1, Pt3, Pt4, Plane);
  166.     printf("[%lf %lf %lf %lf %lf]\n",
  167.        Plane[0], Plane[1], Plane[2], Plane[3], Plane[4]);
  168.  
  169.     TrivPlaneFrom4Points(Pt1, Pt1, Pt1, Pt4, Plane);
  170.     printf("[%lf %lf %lf %lf %lf]\n",
  171.        Plane[0], Plane[1], Plane[2], Plane[3], Plane[4]);
  172.  
  173.     TrivPlaneFrom4Points(Pt1, Pt1, Pt1, Pt1, Plane);
  174.     printf("[%lf %lf %lf %lf %lf]\n",
  175.        Plane[0], Plane[1], Plane[2], Plane[3], Plane[4]);
  176.  
  177.     for (i = 0; i < 100; i++) {
  178.     for (j = 0; j < 4; j++) {
  179.         Pt1[j] = IritRandom(-10, 10);
  180.         Pt2[j] = IritRandom(-10, 10);
  181.         Pt3[j] = IritRandom(-10, 10);
  182.         Pt4[j] = IritRandom(-10, 10);
  183.     }
  184.  
  185.     TrivPlaneFrom4Points(Pt1, Pt2, Pt3, Pt4, Plane);
  186.     if (!APX_EQ(-Plane[4], (Plane[0] * Pt1[0] + Plane[1] * Pt1[1] +
  187.                     Plane[2] * Pt1[2] + Plane[3] * Pt1[3])) ||
  188.         !APX_EQ(-Plane[4], (Plane[0] * Pt2[0] + Plane[1] * Pt2[1] +
  189.                     Plane[2] * Pt2[2] + Plane[3] * Pt2[3])) ||
  190.         !APX_EQ(-Plane[4], (Plane[0] * Pt3[0] + Plane[1] * Pt3[1] +
  191.                     Plane[2] * Pt3[2] + Plane[3] * Pt3[3])) ||
  192.         !APX_EQ(-Plane[4], (Plane[0] * Pt4[0] + Plane[1] * Pt4[1] +
  193.                     Plane[2] * Pt4[2] + Plane[3] * Pt4[3])))
  194.         printf("Error (%d) [%lf %lf %lf %lf %lf]\n", i,
  195.            Plane[0], Plane[1], Plane[2], Plane[3], Plane[4]);
  196.     }
  197. }
  198.  
  199. #endif /* TEST_PLANE_4D */
  200.  
  201.